Generalized Pomeranchuk instabilities in graphene 
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We study the presence of Pomeranchuk instabilities induced by interactions on a Fermi liquid 
description of a graphene layer. Using a recently developed generalization of Pomeranchuk method 
we present a phase diagram in the space of fillings versus on-site and nearest neighbors interactions. 
Interestingly, we find that for both interactions being repulsive an instability region exists near 
the Van Hove filling, in agreement with earlier theoretical work. In contrast, near half filling, the 
Fermi liquid behavior appears to be stable, in agreement with theoretical results and experimental 
findings using ARPES. The method allows for a description of the complete phase diagram for 
arbitrary filling. 
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INTRODUCTION 

Correlated electron systems in two dimensions have at- 
tracted a lot of attention in the last years, especially 
due to an important number of experiments that pro- 
vide undisputable evidence of the existence of new exotic 
phases of matter. 

One such example corresponds to nematic and stripe 
(smectic) phases in high T c superconductors in the under- 
doped region and fractional quantum Hall effect systems 
at high magnetic fields [lj. A nematic phase is char- 
acterized by orientational but not positional order and 
it has been proposed to explain the observed transport 
anisotropics. One important point about these phases is 
that they arise spontaneously, decreasing the rotational 
symmetry without a lowering of the lattice symmetry. 
Another more recent case is given by strontium ruthcn- 
ate Sr3Ru207, which is well modeled as a bilayer system 
and shows a large magnetoresistive anisotropy. This ob- 
servation has been argued to be consistent with an elec- 
tronic nematic fluid phase. Experimentally, two consec- 
utive metamagnetic transitions have been observed and 
the region in between has been proposed to be a conse- 
quence of a Pomeranchuk instability, due to a nematic 
deformation of the Fermi surface, in very close analogy 
to what happens in fractional quantum Hall gallium ar- 
senide systems 0]. Yet another interesting material is 
the heavy fermion compound URu2Si2 in which a hidden 
order phase arises through a second order transition at 
around 17.5°K. The order parameter of this new phase 
has remained elusive to theorists up to date. Different 
types of order have been proposed, but the situation is 
still under debate In recent work, Varma and Zhu 
Q have proposed that this phase transition could corre- 
spond to a Pomeranchuk instability inducing a deforma- 
tion in the antisymmetric spin channel, stabilized by a 



phase characterized by a helicity order parameter. 

The experimental findings mentioned above triggered 
different theoretical studies on low dimensional corre- 
lated systems to search for such exotic phases [J. Spe- 
cial attention has been paid to the possibility of Pomer- 
anchuk instabilities [51 givi ng rise to such novel phases 

0, SJp, E3, EL Hill, 13, M, Gl- In a previous 

paper [171 ] . motivated by these investigations, we devel- 
oped a generalization of Pomeranchuk's method to search 
for instabilities of a Fermi liquid. The method we pre- 
sented is applicable to any two dimensional lattice model 
with an arbitrary shape of the Fermi surface (FS) at zero 
temperature. The main results of our previous paper 
were summarized in the form of a recipe whose steps we 
give below for completeness. Our method is particularly 
well suited to analyze systems with weak interactions and 
then graphene appears as a perfect arena to test our tech- 
niques, since electron interactions are argued to be small 
due to strong screening. 

The recent isolation of graphene (HI, the first purely 
two-dimensional material, which is made out of carbon 
atoms arranged in a hexagonal structure, led to an enor- 
mous interest and a large amount of activity in study- 
ing its properties. The low doping region, near to half 
filling, became the subject of attention due to the pe- 
culiar behavior described by chiral massless charge car- 
riers. Several Van Hove singularities are present at en- 
ergies of the order of the hopping parameter E ~ 2.7 
eV and these singularities are expected to have an im- 
portant role in the properties of the system. Although 
in first approximation graphene layers are well modeled 
by free fermions hopping on a hexagonal lattice, there 
have been a number of works in the literature where the 
effects of electron-electron interactions were taken into 
account 0,11. However, such analysis were centered 
in the undoped (half- filling) regime or very close to it, 
and explicit analytic results for a wide range of fillings 
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are still lacking. 

In the present paper we apply the method developed 
in [ItJ to fermions in a graphenc sheet in the presence 
of electron-electron interactions up to nearest neighbors. 
A previous work using a mean field approach, showed 
that long range interactions could lead to Pomeranchuk 
instabilities 21[ at the Van Hove filling. On the other 
hand, at dopings very close to half-filling, it was antici- 
pated that graphene should behave as a Fermi liquid [12] • 
This was later confirmed experimentally by ARPES ex- 
ploration of the FS 23[. The results presented here are 
consistent with these findings, while our method allows 
for a more complete and systematic study of the whole 
space of fillings. 



TWO DIMENSIONAL POMERANCHUK 
INSTABILITIES: REVIEW AND 
IMPROVEMENT OF THE METHOD 

In this section we will review the generalization of 
Pomeranchuk method first presented in 17], and dis- 
cuss a shortcut that can be used as an alternative to 
the change of variables proposed there. 



S (k) 



(2) 



Review of the method 

According to Landau's theory of the Fermi liquid, the 
free-energy E as a functional of the change Sn^ in the 
equilibrium distribution function at finite chemical po- 
tential fx can be written, to first order in the interaction, 
as 



E = Jd\(e(k)-tj,)5nk+ ~ Jd\ 7dV/(k, k') 



SnkSnu 



(1) 



Here e(k) is the dispersion relation that controls the free 
dynamics of the system, the interaction function /(k, k') 
can be related to the low energy limit of the two particle 
vertex. Note that we are omitting spin indices, and con- 
sidering only variations of the total number of particles 
?ik = ?T,|k + fij.k- This implies that the considerations 
that follow will be valid in the absence of any external 
magnetic field and at constant total magnetization. The 
identification of these two functions is the starting step 
of our calculation: 

Step 1: write the energy as in |T]) and identify the func- 
tions e(k) and /(k, k'). 



Instead of the cartesian variables in momentum space 
(k x ,k y ) we will find convenient to define a new set of 
curvilinear variables (<?, s) according to 

.9 = ff(k) = M-e( k ) , 



The variable g varies in the direction normal to the Fermi 
Surface (FS), whose position is defined by g = 0. We 
choose s such that it is constant at constant distance to 
the FS, varying in the longitudinal direction tangent to 
the FS, namely it satisfies the restriction 



Vs(k).Vg(k) = 0. 



(3) 



Since solving this for s may be a difficult task, we de- 
velop bellow an alternative procedure that replaces this 
calculation. Even if not mandatory, we will chose the 
variable s such that we give a complete turn around each 
connected piece of the FS when it runs from — n to n. 

From the above change of variables we obtain the Ja- 
cobian 



d(g, s) 



(4) 



which is the relevant outcome of this step of the calcula- 
tion: 



Step 2: with the help of the dispersion relation e(k) 
obtained in step 1, change the variables according to (|2|) 
to obtain the Jacobian 

In a stable system, the energy ([1]) should be positive for 
all (5?ik that satisfy the constraint imposed by Luttinger 
theorem 24 1 



0. 



(5) 



Pomeranchuk's method roughly consists on exploring 
the space of solutions of constraint j5]) to find a Sn^ that 
turns the energy into negative values, thus pointing to 
an instability of the system. 

In terms of our new variables g and s, we can write 
<5n k ( g s ) at zero temperature as 



= H[g + Sg(g,s)]-H[g] 



(6) 



where H is the unit step function and 8g(g, s) is an small 
perturbation parameterizing the deformation of the FS. 
Replacing into the constraint ([5]), changing the variables 
of integration according to ([2]) and performing the inte- 
gral to lowest order in 5g we get 



dsJ(s)6g(s) = . 



(7) 



Here J(s) = J(g, s)\ g= o and Sg(s) = 5g{g, s)| 5= o- In case 
the FS has a nontrivial topology, the integral includes a 
sum over all different connected pieces. 

We see that in order to solve the constraint 5g(s) can 
be written as 



<fo00~ J- 1 (s)d s \(s) , 



(8) 
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in terms of a free slowly varying function A(s). Even 
if in (|5|) and ([7]) a sum over different connected pieces 
of the FS may be assumed, this particular solution does 
not consider excitations in which some particles jump 
between different connected pieces . 

Using the change of variables ^ and with the help of 
eqs. © and ([5]), we can rewrite the energy E to lowest 
order in Sg as 

E = J ds'J ds^(s')\ (j-\s)8(s - s 1 ) + f(s,s')\ 4(a) , 

(9) 

where we call d s X(s) = V'( s ) an d f(s,s') = 
f(9,s;g',s')\g=g, =0 . 

Note that the left hand side of the stability condition 
E > has two terms, the first of which contains the in- 
formation about the form of the FS via J _1 (s), while 
the second encodes the specific form of the interaction in 
f(s,s'). There is a clear competition between the inter- 
action function in the second term of the integrand and 
the first term that only depends on the geometry of the 
unperturbed FS. 

We see in ([9]) that the energy E is a bilinear form 
acting on the real functions ip(s) that parameterize the 
deformations of the FS 

E=(1>,4). (10) 

This is the next step of our calculation that can be sum- 
marized as: 

Step 3: write the energy as the bilinear form j9]) using 
the functions J(s) and f(s, s) identified in steps 1 and 2. 

The stability condition is then equivalent to demanding 
this bilinear form to be positive definite for any possible 
smooth deformation of the FS 

W>€ L 2 [FS] : (V>,^>) > 0, (11) 

where L2PS] stands for the space of square-integrable 
functions defined on the Fermi surface. 

In consequence, a straightforward way to diagnose an 
instability is to diagonalizc this bilinear form looking for 
negative eigenvalues. To see that, we choose an arbitrary 
basis of functions {ipn}neN of L2PS], in terms of which 
we can write 

^00 =$>„Vn(s) , (12) 

n 

and the the stability condition now reads 

E = ^ O-nami^njtpm) > . (13) 
n, m 



This configures our 

Step 4: choose an arbitrary basis {ip n }n£N of the space 
of functions L 2 [FS] . 

The above defined bilinear form can be considered as a 
pseudo-scalar product in L 2 [FS] . In general the functions 
of the basis {tpn}neN will not be mutually orthogonal 
with respect to this product. Moreover in the presence 
of instabilities, the pseudo-scalar product may lead to 
negative pseudo-norms (V','0) < 0. 

We can then make use of the Gram-Schmidt orthogo- 
nalization procedure to obtain a new basis of mutually 
orthogonal functions {(„}„ e n . In terms of them an arbi- 
trary deformation of the FS parameterized by a function 
■0(s) can be decomposed as 

#0=E b nUs), (14) 
n 

which implies that the stability condition on such defor- 
mation will read 

# = E im 2 >0 - ( 15 ) 

n 

In summary: 

Step 5: apply the Gram-Schmidt orthogonalization pro- 
cedure to go from the arbitrary basis {ipn}neN chosen 
on step 3 onto a basis of mutually orthogonal functions 

In (|15|) we note that the only possible source of a neg- 
ative sign is in the pseudo-norms Xn = In case 
the i-th pseudo-norm \i is negative, a deformation pa- 
rameterized by ip{s) tx £i(s) is unstable. In other words 
the pseudo-norms {Xn}neN can be taken as the stability 
parameters, a negative value of Xi indicating a instability 
in the i-th channel. Then: 

Step 6: compute the pseudo-norms of the new basis 
functions {£n}n6 n- If for a given channel i one finds that 
Xi = (&>&) < 0, the FS is diagnosed to be unstable. 

These six steps summarize the generalized Pomer- 
anchuk method. It can be applied to any two dimensional 
model with arbitrary dispersion relation and interaction. 
Note that since L 2 [FS] is infinite dimensional, the present 
method is not efficient to verify stability: at any step i it 
may always be the case that for some j, Xi+j < 0. More- 
over, in the case of nontrivial topology, excitations con- 
sisting on particles jumping between different connected 
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pieces of the FS were not included in the solution of the 
constraint ([7]) , and they may lead to additional instabili- 
ties. The same is true excitations involving spin or color 
flips, that we arc not considering. 



An alternative to the change of variables 



that in order to be compatible with ([18)1 imply 

|V ff (k( s ))| 



J~ L (s) = 



|k(.s)| 



(24) 



Then with all the above, wc can replace our previous 
step 2 by a new version 



As advanced, to avoid the task of solving the constraint 
((H]) that defines the variable s, we will develop here an 
alternative procedure to derive the form of the Jacobian 
evaluated on the FS J(s). 

We begin by defining a paramctrization of the FS 



k(t) = (k x (t),ky(t)) , 



-7T<<<7T, (16) 



in terms of an arbitrary parameter t. In other words, 
given the function g(k) defined in ([2]), we choose k(t) such 
that Vi : g(k(t)) = 0. In terms of such paramctrization 
we can decompose the Dirac 6 function as 



%(k)) = / dt 



|k(i)l 



|Vff(k(t))| 



k-k(i)), (17) 



(19) 



(a proof of this formula is given in the Appendix) . 

The integral I of an arbitrary function F(g,s) along 
the FS can be written as 

/ = J dsF(g lS )\ g=0 , (18) 

or in other words 

I = j dsdgF(g,s)S(g). 

Changing variables to k 

1 = J d 2 kF( S (k), 5 (k)) J- 1 ( S (k), 5 (k))%(k)), (20) 

and replacing p7|) we get 

/ = J <fkF(s(K),g(k)) J-\s(k),g(k)) x 

X I dt \M))/ 2){k - m) ' (21) 

or, interchanging the order of the integration and per- 
forming the k integrals 

I = JdtF(s(k(t)),g(k(t))) J- 1 ( 8 (k(t)),g(k(t))) lv l gHf m ■ 

(22) 

Now using the fact the g(k(t)) = and defining the pa- 
rameter t such that s(k(t)) = t, we get 



ds F(s,g)\ J- 1 (a) 



|k( S ) 



|V.g(k( S ))| 



(23) 



Step 2': with the help of the dispersion relation e(k) 
obtained in step 1, parameterize the FS and obtain the 
Jacobian evaluated on the FS according to ([24]) . 



Then, even if it may be very difficult to solve the par- 
tial differential equation (J3j) in order to explicitly obtain 
the Jacobian, its restriction to the FS is all what we need, 
and can be obtained by the much easier task of parame- 
terizing the FS. 



POMERANCHUK INSTABILITY IN GRAPHENE 

In the present section wc will apply the method re- 
viewed in the previous pages to the case of fcrmions in a 
graphenc layer with Coulomb interactions. 

Free Hamiltonian: tight-binding approach 

Graphene is made out of carbon atoms arranged in a 
hexagonal lattice. It is not a Bravais lattice but can be 
seen as a triangular lattice with two atoms per unit cell. 
The tight-binding Hamiltonian for electrons in graphenc 
considering that electrons can hop only to nearest neigh- 
bor atoms has the form 



H = -t ( a lAi + H - c 



(25) 



where a a i , b a i are the creation and annihilation operators 
related to each of the unit cell atoms. 

The diagonalized Hamiltonian can be written in terms 
the occupation numbers of rotated lattice operators de- 
fined by 



.k^Kki^^k) (& CT k±9^a CTk ) , (26) 



h*(k) 



,t 



e(k) 



h(k) 



e(k) 



where the function h(k) satisfies |/i(k)| 2 = e(k) 2 and 
reads 



h(k) = t \ cos{k x ) - 2ism.(k x -l) + 2cos 2 (y) + 

+4i cos(— p^)sin(y) +4cos(— ^-)cos(fc x ) 



(27) 
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It results in 

H = f d 2 kJ2 (4(k)n+ k + s-Qt)n- k ) , (28) 

here we defined the energy bands by (k) = ±e(k) with 

e(k) = t^j 3+4 cos ^fc^a^ cos ^^-k y a^j + 2 cos (V3k y a) , 

(29) 

where a is the carbon-carbon distance (a ~ 1.24 A). 

The energy resulting from a small variation on the oc- 
cupation numbers on ([28)1 at finite chemical potential fi 
reads 



E = J A( £ (k)-/i)(<5 



"k - <K ) i 



(30) 



where use have been made of the fact that the dispersion 
relation do not have spin indices to write the resulting 
expression in terms of variations of n k = + n, k . 

To fix the ideas we consider only non- vanishing varia- 
tions of the occupation numbers in the upper band, i.e. 



Sn,, 



0, 5rq ^ 0. This can be done without loss of 



generality due to the symmetry of the system under the 
interchange of the positive and negative bands. We then 
have 



E Q = Jd 2 k (e(k) - n)8n+ . 



(31) 



This has the form o the free term of ([T]), giving us one 
of the ingredients required by our step 1 defined above, 
namely the dispersion relation e(k). 

The band structure is shown in Fig[TJi. The density of 
states is presented in Fig. [TJd where a Van Hove singular- 
ity can be seen to be present for /i = ±t. The resulting 
FS's for different fillings are shown in FigQJ:. 



Interactions in graphene 

To complete the ingredients required by our step 1, we 
need the quasiparticle interaction function /(k, k'). In 
what follows, for completeness and to set up our conven- 
tions, we briefly describe how to derive its exp ression to 



first order in a pcrturbative expansion 

We will consider density-density interactions, both on- 
site (with strength U) and between nearest-neighbors 
(with strength V), namely our interaction Hamiltonian 
reads 



U 



Hi n t = — ^V t (n. ( - 1) 



V niHj 



(32) 



where stands for nearest neighbors, and the den- 

sity operators rij = n,j + refer to the original lattice 
operators , bi . 




A 














V/ 



FIG. 1: (Color online) a) Energy spectrum for the tight- 
binding approach, b) Density of states per unit cell as a 
function of the energy. All the quantities are given in units of 
t. c) Left: FS for t < fj, < 3i. Right: FS for < fx < t 



One then computes the energy in the mean field ap- 
proximation. The result can be written in terms of the 
mean field values of the occupation numbers of the ro- 
tated lattice operators that diagonalize Hq, reading 

k,k',a,/3 

x (Uaip + j (F(0)(<x° 5 +<7^))) - 

lu J]] ( n ka n k'/3+ n kQ n k'/3 — n ka n k'/3 — n ko, n k' 1 fl) X 



2N 

k,k',a,/3 

x^^F(k-k)^Mk') 



h*(k)h(k') J ' 
where N is the number of sites and 



(33) 



ik-<5„ 



with 



a=l 



*i = o(i,^), 
& = a(±,-^) 
«a = o(-l,0). 



(34) 

As advanced in the previous section, we will concen- 
trate in variations of the occupation numbers that keep 
the total magnetization constant. In other words, we 
assume Sn^ k = 5n^ k . Similarly to the free part, the in- 
teractions between quasiparticles can then be written in 
terms of the variation of the total number n k = n-^ k +n^ k 
as 

E int = Jdkdk! f(k 1 k')6n+6n+ . (35) 
The function /(k, k') is then obtained from the mean 
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field value of the energy 



(36) 



G{x) = \(j2 - 2cos(V3x) - sec^- 



Low energy sector: < t 



V3c 



(39) 



We have then completed step 1, getting the dispersion 
relation (|29|) and the interaction function (f36|) . 



Parametrization of the Fermi surface 

Step 2' requires the parametrization of the FS, for 
which we need to study separately fillings lying above and 
below the Van Hove filling. In what follows we present 
the parameterized curves used throughout the calcula- 
tions. 

High energy sector: > t 

We call high-energy sector the case in which t < < 
3t, i.e. the region of fillings which lie above the Van Hove 
singularity. As can be seen in Figfl}:, for fi/t ~ 3 the FS 
is approximately circular, while for values closer to 1 the 
FS takes a hexagonal form. 

In this sector the FS can be parameterized as follows 



k ff ( S ) = (fcf( S ), fcf( S )) , 
where 



-7T < s < 7T , (37) 



00 = ^sign(s) arccos[G(fcf (s))} 



3 a 
1 



\/3a 



i>H - uh \s\) ■ 



(38) 



with loh, 4>h and the auxiliary function G(x) defined as 



luh = — arccosl 

7T 



2 arccos , 



The low energy sector corresponds to fillings satisfying 
< \fi\ < t. In this case the FS consists of pockets 
centered at the vertices of the Brillouin zone, as can be 
seen in Fig[T}:. By using the periodic identifications of the 
momentum plane, we see that only two of them are non- 
equivalent. In consequence one can describe the total FS 
as two pockets centered in the two Dirac points k± = 

(o,±&). 

For example, for the FS pocket centered in (0, we 
can choose a parametrization of the form 



k L ( S ) = (k^( S ), k$ (8)) , 



-IT < S < 7T , (40) 



with 



fc£( s ) 



— sign(s) arccos [G(fc^(s))] 
i 

-ul N), 



\/3a 



(41) 



and 



—t—fj, 



— f arccosl 

7T V \ 2t 



ffi-t 

V 2t 



2 arccosl ■ 



-t—fj, 



V 2t 



(42) 



With the above parametrizations of the high and low 
energy sectors, we can compute the corresponding Jaco- 
bian evaluated on the FS according to (|2"4")) . obtaining 



J-\ S ) 



^g nv /6 M 2 - Ai 4 + 4(^_i) cos (^(f-| s |))-2cos(^(7r-2| S |))-3, t< | M | <3i, 



^g^6/i 2 - A t 4 + 4 JJ? - 1) cos (0l-wl|s|) - 2 cos (2 {^ l -lo l \s\)) - 3 , < |m! < * ■ 



(43) 



This completes our step 2', providing us with the values 
of the Jacobian evaluated on the FS J~ 1 (s). The interac- 
tion function evaluated on the FS f(s, s') is obtained by 



simply replacing the paramcterizations of the high and 
low energy sectors in (f3"6")) . This allows us to complete 
step 3, by constructing the energy function as a bilinear 
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U 



-10 



V 




FIG. 2: (Color online) 3D Phase diagram for graphene. The 
phase diagram was constructed by exploring the first 20 
modes, the shaded region is unstable. The first three modes 
cover most of the instability region, the remaining modes just 
re-draw the details of the boundary. Note that for purely re- 
pulsive interaction (positive U and V) there is an unstable 
region near the Van Hove filling fi ~ I . On the other hand, 
near half fillings fi ~ the Fermi liquid is stable for any value 
of the interactions. 



form. 



Instabilities and phase diagram 

To proceed with step 4, we choose a basis of the 
space of functions L2 [FS] . Here for simplicity we choose 
trigonometric functions 

{-0n(s)}„ e w = {cos(tos), sin(ms)} meA r ■ (44) 

Following step 5 by means of an orthogonalization pro- 
cedure, we obtain the new basis of mutually orthogonal 
functions {£„}„ G jv- 

In terms of this new basis and according to step 6, we 
compute the stability parameters Xi = (£»)&}• These are 
functions on the space of parameters (/i, U, V) that may 
become negative in some regions. If this is the case, in 
such regions the Fermi liquid is diagnosed to be unstable. 

In our calculations, this last step was performed nu- 
merically, due to the complication of the integrals in- 
volved in the pseudo-norms {Xn}neN- 



RESULTS 

The phase diagram on the space spanned by the in- 
teraction strengths U, V and the chemical potential fi 
is shown in Fig [5J There we plot the instability region 
determined by the dominant unstable channels. 

The method used in this work allows to explore all pos- 
sible fillings and to draw a phase diagram for graphene 
systems valid both below and above half-filling. The ad- 
vantage of our approach lies in the fact that it can be 
pursued systematically following the steps described in 
previous sections, studying separately each deformation 
mode of the FS. 

For fillings around the Van Hove filling the Pomer- 
anchuk instability is favored. For on-site and nearest- 
neighbor Coulomb repulsion (U > and V > 0) we find 
a region of parameter space where the system presents an 
instability. Near to the Van Hove filling our results are in 
agreement with those found using a mean field approach 
in (23. 

On the other hand, we do not see any instability 
around half filling. This is in agreement with the re- 
sults presented by Sarma et al in [22| for doped graphene, 
where the authors found that extrinsic graphene is a well 
defined Fermi liquid for low energies, within the Dirac 
fermion approximation. Moreover, this agrees with ex- 
perimental results found using ARPES presented in [23[ , 
implying that graphene is a Fermi liquid for low dopings. 

For attractive on-site interaction, the region where the 
instability is detected depends strongly of the nearest 
neighbors interaction strength. 

We find that even the smoothest deformations of the 
FS, i.e. those described by lower modes in our orthogo- 
nal basis, present instabilities. Indeed, they cover most 
of the unstable region. In Fig[3] instability regions cor- 
responding to the first modes are drawn. In Fig[3K the 
V = plane is shown, the unstable mode corresponding 
to the colored region is the 0-th mode and the correspond- 
ing FS deformation is presented in Fig[3J:. In Fig[3)3 and 
[3ji the instability regions for the first modes are plotted 
in the U-V and fi-V plane respectively. In both Fig- 
ures the xo instability can be seen, and a new region 
appears where an instability of the \i and X2 channels 
is present. This instability appears for positive values of 
the interaction strength and is closer to the Van Hove 
filling. This is in agreement with earlier results, where 
a Pomeranchuk instability at the Van Hove filling was 
found within a mean field approximation [2lj . The inter- 
acting FS presented there has the same geometry than 
that of the corresponding deformation channel \2 shown 
in the left column of Figj3j:. 

The results presented above confirm earlier predictions 
about the Fermi liquid behavior of graphene for the cases 
near Van Hove filling and near half filling. They also pro- 
vide a more complete description of the phase diagram 
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FIG. 3: (Color online) Instability sectors for the first unstable 
modes, a) Instability for the xo mode in the V = plane. It 
is reached only for attractive on-site interactions (U < 0). b) 
Instability regions on the plane fi = 0.99: the modes xi an d 
Xi are unstable for U > and V > 0. c) FS deformations 
corresponding to the first unstable modes: the left column 
corresponds to the low energy sector and right column to the 
high energy one. d) Instability reached by the xo, Xi an d 
X2 modes in the U = plane. Remarkably, it is present for 
repulsive interactions (U > 0,V > 0). 



The efficiency of the method shows up in the fact that 
the complete phase diagram is obtained by exploring a 
few number of modes. The introduction of higher modes 
does not enlarge substantially the region of instability 
but re-draw the details of its boundaries. 

The method can be applied either analytically or nu- 
merically according to the complexity of the system un- 
der investigation. In the previously studied case [T3] , the 
analytical calculations were pursued up to the end, al- 
lowing us to draw the phase diagram. In the present 
case, the calculations were performed analytically up to 
the point of the evaluation of the instability parameters 
Xi, which involved complicated integrals that were then 
performed numerically. The method is suitable for di- 
rect application to numerical data encoding the form of 
the Fermi surfaces, like those obtained by application of 
ARPES. 

The form of the method presented here is suitable for 
any two-dimensional lattice model at zero temperature. 
However, it does not consider instabilities arising from 
particles jumping between different disconnected pieces 
of the FS or from spin or color flipping. It can be easily 
extended to consider these effects, as well as to three di- 
mensional systems, such as multilayer graphene, ruthen- 
ates, etc. The generalization to finite temperatures in- 
volves a different definition of the pseudo scalar product. 
All these extensions will be presented in a forthcoming 
work 
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for the entire range of fillings in between these two lim- 
iting cases. 



SUMMARY AND OUTLOOK 

We have explored Pomcranchuk instabilities in 
graphene using a recently developed generalization of 
Pomcranchuk method. We have obtained the phase di- 
agram of the dominant instability as a function of the 
on-site U and nearest neighbor V interaction strengths 
and the chemical potential fi (Fig[2]). We analyzed sev- 
eral planes of the 3D phase diagram obtaining a good 
agreement with previous theoretical results and experi- 
mental findings (Fig|3j) . 

The phase diagram makes apparent some interesting 
features of the system. For example no instability is de- 
tected at low energies. This behavior is noteworthy be- 
cause in this sector a Dirac massless fermions approach 
can be used to describe the graphene layer in the absence 
of interactions. On the other hand, for energies close to 
the Van Hove energy, instabilities appear to cover a large 
region in the U-V plane. 



Given the implicit curve defined by .g(k) = 0, we can 
choose a parametrization in terms of a vector function 



k(t) = (k x (t),ky(t)) , 



-7T < t < TT . 



(45) 



depending on an arbitrary parameter t. and defined so 
that Vt : <?(k(t)) = 0. In terms of this parametrization 
we want to prove that the Dirac S function can be written 
as 

^ = /*jvS&^- k ^' (46) 

To that end, we write more explicitly the right hand side 



%(k)) = / dt 



|Vff(k(t))| 



6(k x -k x (t))6(k y -k y (t)), 

(47) 

and then rewrite the k x delta function using the well 
known one dimensional formula 

S(f(x)) = 6{ l;_f where f{x) = , (48) 



/'(*) 



to get 



'(900) = / dtJ^L ^T^ 11 S(k y - k y (t)) , 



'|V fl (k(t))| k x (t) 



(49) 
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where we call t{k x ) to the solution of k x — k x (t) = 0. 
Performing the t integral 



<%(k)) = 



|k(fc x )| 5(ky~k y (k x )) 



|Vff(**)l k x (k x ) 



(50) 



where we use the notation f{k x ) = f{t(k x )) for any func- 
tion /. Writing explicitly the square roots on the vector 
norms 



<%(k)) 



■sj ' h%{hx) -\- h^hx) S{ky ky (k x )) 
y/ (d x g(k x , k y (k x ))) 2 + (d y g(k Xl k y (k x ))) 2 k x (k x 



(51) 



A further rearrangement of the formulas gives 
%(k)) = 



V^+(^) a *(*»-*»(*»)) 



1 



a yS (fc x ,fc H (fc x )) ) °v9[ K x, K y {K x )) 



(52) 

where we can identify the derivative in the numerator 
with the quotient of derivatives in the denominator to 
cancel the square roots, obtaining 



%(k)) 



5{ky ky{kx)^ 



(53) 



dyg{k Xl kyikj^)^ 

if k y takes the place 



This is an identity in virtue of 
of x. 

Then we just proved formula (|47p that we used during 
our calculation of the Jacobian of the change of variables 
evaluated in the FS. 
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